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Abstract 

We consider the linear Schrodinger equation and its discretization by split-step 
methods where the part corresponding to the Laplace operator is approximated 
by the midpoint rule. We show that the numerical solution coincides with the 
exact solution of a modified partial differential equation at each time step. This 
shows the existence of a modified energy preserved by the numerical scheme. 
This energy is close to the exact energy if the numerical solution is smooth. As a 
consequence, we give uniform regularity estimates for the numerical solution over 
arbitrary long time. 
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1 Introduction 

We consider the linear Schrodinger equation 

dtu(t, x) = — iAu(t, x) + iV(x)u(t, x), u(0, x) = u°(x), (1.1) 

with initial condition and potential function V(x) G R. The wave function 
u(x, t) depends on x € T d or M. d and the time t > 0. The operator A is the 
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d-dimensional Laplace operator. In the following, we consider mainly the case 
where x S T d . The case of the whole space is totally similar. The equation (1.1) 
is symplectic and its solution preserves the I? norm and the energy 

u>-> \Vu\ 2 + V\u\ 2 dx = (u\ - A + V\u). (1.2) 

The solution of (1.1) is given by 

u(t, x) = exp(it(-A + V))u°(x), 

and a standard method to simulate this solution is to consider the approximation 

exp(ih(-A + V)) ~ exp(-i/iA) exp(ihV) (1.3) 

for a small stepsize h > 0. The solution at a given time t = nh is then approxi- 
mated by 

exp(ii(— A + V))u ~ I exp(— ihA)exp(ihV)Yu°. (1.4) 

The advantage of this method is that it yields a symplectic scheme preserving 
the L 2 norm. Moreover, it is very easy to implement by using the fast Fourier 
transform: while the operator A is diagonal in the Fourier space, the operator V 
acts as a multiplication operator in the phase space. For finite time, this splitting 
scheme yields a consistent numerical scheme: as h — ► and if the numerical 
solution is smooth, it can be shown that (1.4) yields a convergent approximation 
of order 1 in h, see [12]. Considering higher order approximation such as the 
symmetric Strang splitting or higher order splitting methods allows to obtain 
higher order approximation scheme under the assumption that the numerical 
solution is smooth enough, see [12, 9]. 

Concerning the long-time behaviour of such methods, very few results exist. In 
[3], Dujardin & Faou showed the conservation of the regularity of the numerical 
solution (1.4) in T 1 over very long time, provided the potential function is small 
and smooth. Moreover, even in this situation, resonances effects appear for some 
values of h: typically when exp(— ihA) posseses eigenvalues close to 1. 

In the finite dimensional case, the long time behaviour of splitting method 
can be understood upon using the Baker-Campbell-Hausdorff formula (see for 
instance [8]). Roughly speaking, this result states that for two matrices A and 
B, we can write 

exp(L4)exp(tB) = exp(tZ(t)) 

where Z{t) = A + B + t[A, B] + t 2 ■ ■ ■ , with [A, B] = AB - BA the matrix com- 
mutator. Hence the long time behaviour of the numerical solution corresponding 
to (1.4) can be analyzed by considering the properties of the matrix Z(t) which 
is a small perturbation of the original operator A + B for small time t. However, 
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to be valid, the BCH formula requires h to be small enough with respect to the 
inverse of the norms of A and B. This makes this strategy impossible to apply 
directly for unbounded operators, unless a drastic CFL like condition is used for 
the full discretization of (1.1). 

In this paper, we consider the time discretization 

exp(i/i(-A + V)) ~ exp(ihV)R(-ihA) (1.5) 

where 

l + z/2 

R{Z) = 

is the stability function of the midpoint rule. Such an approximation is clearly 
consistent with (1.1) if the solution is smooth enough. Moreover, it defines a 
symplectic numerical scheme preserving the L 2 norm, and easily implemented 
by using the fast Fourier transform. Similar schemes have been considered in 
[1, 13, 16]. 

Recall that for all x S R we have 

— %X = exp(2i arctan(x)). 
1 — ix 

and hence we can write 

, , 1 - ihA/2 . , hA n 

R(-ihA) = - - . hA ^ = exp(2*arctan ( - — )), 

where now 2arctan ( — A^) is a bounded operator from L 2 to itself. Using this 
representation, we show in this work that there exists a symmetric operator S(h) : 
L 2 — > L 2 such that 

ex.p(ihV)R(-ihA) = exp(ihS(h)), 

with 

S(h) = -j- arctan {^—) + V(h) 

where V(h) : L 2 —> L 2 is a modified potential. 

Hence, for all n and all initial value u°, we have 

u n = (exp(ihV)R(-ihA)) n u° = exp{inhS{h))u° 

and hence the numerical solution u 11 coincides with the exact solution of the 
modified equation 

dtu = S(h)u 

at each time step t n = nh. This implies that the associated energy 

(u | S(h) | u) 
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is preserved along the numerical solution associated with the split-step scheme 
(1.5). Moreover this energy is close to the original energy (1.2) if u is smooth. 
Using these properties, we give regularity bounds for the numerical solution over 
arbitrary long time. 

Such a result is to our knowledge the first extension in an infinite dimen- 
sional setting of the classical backward error analysis for Hamiltonian ordinary 
differential equation (see [8, 11]). Note in particular that as in the case of linear 
ordinary differential equation, this result is valid for arbitrary long time, while 
such results classically hold for exponentially long time with respect to the step 
size for nonlinear ordinary differential equations. 

It is worth noticing that such result does not hold hold for the splitting scheme 
(1.3) for which it is known that resonance effects occur, see [3]. The main differ- 
ence between (1.5) and (1.3) lies in the high frequencies regularization effect of 
the midpoint rule: by essence, the logarithm of the operator R(—ihA) is bounded 
while the logarithm of exp(— ihA) is not well defined when hA possesses eigen- 
values close to multiples of 2ir. Note that this does not affect the approximation 
property of the scheme for finite time and smooth numerical solution. 

Similarly this result does not automatically extend to situations where the 
propagator R(—ihA) is replaced by a higher order approximation of exp(— ihA), 
or for higher order splitting schemes (see [8, Chap III]). We discuss this point in 
the last section of this work, and show by numerical experiments that in general 
resonance effects appear. 

Let us mention that in the nonlinear situation, results exist concerning the 
long-time behaviour of splitting scheme applied to the nonlinear Schrodinger 
equation: see the recent works of Faou, Grebert & Paturel [4, 5] and Gauck- 
ler &; Lubich [6, 7] for the long time behaviour of splitting schemes applied to 
NLS when the initial solution is small. However, to our knowledge no existence 
results for a global modified energy have been proved. Note that in this direction, 
concerning the numerical approximation of solitary wave, Duran & Sanz-Serna 
[2] have proved the existence of a modified solitary wave over finite time for the 
numerical solution associated with the midpoint rule. 



2 Statement of the results 



We represent a function u £ L 2 (T d ) by its Fourier coefficients u 
defined as 
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where for k = (ki,...,kd) G Z d and x = (xi, ■ ■ ■ , Xd) G T d we set k • x = 
fcixi + • • • fcdXrf. We define 

\wf = Yl i nfc i 2 ' and ihi^ = 1^ ( x + \ k \ 2 y\ u k? 

the L 2 and the fP Sobolev norms on T d , where for k = (fci, . . . , G Z d , we set 

For an operator A = (Aj e £) kieZ d acting in the Fourier space C z<i and for a > 1 
we set 

||A|| a = sup|A w |(l + |fc-*| a ). 

We denote by 

C a = {A = {A M ) k/&d | \\A\\ a < oo}. 

If A G £ a with a > d, we can easily show that A G C{L 2 ): see Lemma 4.2 below. 

We say that A is symmetric if for all k,£ G Z d , we have .A^ = or 
equivalently A* = A. In this situation, for w £ L 2 , we set 

(u\ A \ u) = UkAk£U£ = (u, Au) G R 

where ( • , • ) is the L 2 product in T d . For two operators A and I?, we set 

ad A (B) = AS - 5 A 

Finally, with a real function W(x) we associate the operator W = (Wke)k£ez d 
with components W^i = W^-e where W n denote the Fourier coefficient of W 
associated with n G Z d . Thus the operator (Wke)keez d acting in the Fourier 
space corresponds to the multiplication by W. Note moreover that with this 
identification, IIWII < oo with a > d implies that IIWII r _ < oo. 
The goal of this paper is to prove the following results: 

Theorem 2.1 Let a > d, and assume that \\V\\ < oo. There exist hn > and 

7 II II Q, 

a constant C such that for all h G (0, ho), there exists a symmetric operator S(h) 
such that 

exp(ihV)R(-ihA) = exp(ihS(h)), 

satisfying for all h, 

S(h) = - ~ arctan (^-) + V(h) + hW(h) 
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where V(h) and W(h) satisfy 

ii^)ii a +n^)ii a <cmi a , (2-i) 

and where moreover V(h) is given by the convergent series in C a 

V(h) = (dexp Z()W y\v) = V + ^ |f * k *d k Zo{h) (V) (2.2) 

k>l 

with Zo(h) = — 2arctan (—^-), and where the are the Bernouilli numbers. 

Remark 2.2 The size of ho is only proportional to the inverse of ||V"|| , and 
hence is a reasonably small parameter. In particular it does not depend on a 
possible space discretization of the problem through a CFL condition. 

The following result shows that S(h) defines a "modified" energy when applied 
to smooth functions: 

Proposition 2.3 Let (3 G [0,1]. Assume that u G H 1+ ^(T d ), then we have for 
h G (0,ho), 

\(u\S(h)\u) - (u\ -A + V\u)\ < ChP\\uf Hl+fj . (2.3) 
where C depends on [3 and V . 

The next results shows the conservation the modified energy S(h) along the 
numerical solution associated with the split-step propagator. As a consequence, 
we give a regularity bound for the numerical solution over arbitrary long time. 

Corollary 2.4 Assume that u° G L 2 (T d ) and h G (0,h ). For all n > 1, we 

define 

u n = (exp(ihV)R(-ihA)) n u°. 

Then for all n we have 

{u n \S{h)\u n ) = (u°\S(h)\u°). (2.4) 

// moreover u° G H , then there exists a constant Co depending on V and a such 
that for all n G N, 

E l fc | 2 K| 2 + ^ E \<?<C Q \\u%,. (2.5) 

\k\<l/Vh. \k\>l/Vh 

This last result shows that H 1 estimate are preserved over arbitrary long time 
only for "low" modes \k\ < 1/Vh whereas the remaining high frequencies part is 
small in L? . 
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Remark 2.5 The results above obviously remain valid when considering the full 
discretization of (1.1) by collocation methods (see for instance [10]), with esti- 
mates independent of the spectral discretization parameter. 

Remark 2.6 The previous results easily extend to the splitting scheme 

R(-ihA) exp(ihV) 

and to the Strang splitting 

exp(ihV/2)R(-ihA) exp(ihV/2). (2.6) 

Note that in this last situation, the fact that the method is of order 2 allows to take 
[3 £ [0,2] in (2.3). See Section 7 for further details on other possible extensions. 

3 Formal series 

We now start the proof of Theorem 2.1. 
In the following, we set 

Zq := — 2arctan \~^~) 
the diagonal operator with coefficients 

h\k\ 2 

A fc = (Z ) k k = 2arctan {-^-), k G Z d . 

We look for a function t — ► Z(t) taking value into the set of operator acting on 
C zd such that Z(0) = Z and 

Vt€[0,A], e uv e lZ ° =e iZ ®. 

Derivating the equation in t, this yields (see [8]) 

iVe uv e^=i{de WlZ{t) Z'(t))e iZ ®. 

Hence Z(t) has to satisfy the equation (see [8, Chap. III. 4]) 

Z\t) = (dexp^))- 1 ^ = iJ2 ^BA k iZ{t) {V). (3.1) 

fc>0 

and Z(0) = Zq. Here, the are the Bernouilli numbers. Recall that for z € C, 
\z\ < 27r, the expression 

"fcP 



-z 

k>0 



k\ 
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defines a power series of radius 2tt. 
We define the formal series 



z(t) = ^2t i z e 



£>0 

where Zi, I > 1, are unknown operators. 

Plugging this expression into (3.1) we find 

E &- lz t = E t(*E^)V) 

= E^E t* fc E ad^-ad^Cn 

£>0 fc>0 ' €i+-+4=^ 
Identifying the coefficients in the formal series, we find the induction formula: 

V£>1, (£ + l)Z e+1 = ^ jfi h E ad ^i ■■■ ad z ik (V). (3.2) 

fc>0 ' £i+-+£k=£ 

Note that we easily show by induction that for all £, Zg is symmetric. For £ = 1, 
this equation yields 

z i = E f^ ad Un (3-3) 

fc>0 

Note that the main difference with the finite dimensional situation is that the 
"first" term in the expansion is given by an infinite series and that it depends on 
the small parameter h through the operator Zq. The key to control this term is 
to estimate the norm of the operator ad^ . 



4 Proof of Theorem 2.1 

Lemma 4.1 Assume that a > d. There exist a constant C a such that for all 
operator A and B , 

\\AB\\ < C a \\A\\ \\B\\ . 

n n Q "ii iiqH "a 

Proof. We have for k, £ £ Z d , 

\(AB) M \(l + \k- £\ a ) < (l + \k- £\ a ) \A kp \\B kp \ 

l + \k-£\ a 



< \\A 



I \\B\\ y -r- 



(1 + \k -p\ a )(l + \p-£\ a ) 

But as the function x — ► x a is convex for x > 0, we have 

1 + \k - p\ a < 1 + - £\ + \£ - p\) a < 2"" 1 (1 + \k - £\ a + 1 + \£ - p\ a ). 
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Hence we have 

\(AB) M \(l + \k-£\ a ) < 2 Q " 1 || J 4|| ILBII Vf - + ^ ; 

iv Mm "n i;- ii n a n II a Vl+ ib-po 1 + \p-£\ 

and this shows the result, as a > d. ■ 

Lemma 4.2 Let a > d. There exist a constant M a such that for all symmetric 
operator B and for all u G L 2 , we have 

\(u\B\u}\ < M Q ||£|| J|u|| 2 . 

Proof. We have 

\(u\B\u}\ < ^2 \B ke \\u k \\u e \ 



< \\B\\ - — - — -— \uk\\up\ 



< 11*11 E — rr 



k,£ 

2 



using the formula |-Ufc||it^| < ^(|ufc| 2 + |u^| 2 ). This yields the result. ■ 

Lemma 4.3 Recall that Zq = 2arctan an d let W = {Wke)k,iez d be an 

operator. We have for all a > 1 

\\nd Zo W\\ a <<a\\W\\ a . (4.1) 

Proof. For k, i € Z d we have as Zq is diagonal 

(ad Zo W) ke = (Xk-Xe)W ke , 

= (2arctan(/i|A;| 2 /2) - 2 arctan(/i|£| 2 /2)) W u . 
Hence we have for all k, I 6 Z d , 

\( a d Zo W) M \ <ir\W M \ 
and this shows the result. ■ 
Using this Lemma, we see using (3.3) that 

I B I 

ll^ill < WW V -TT-K k < C\\V\\ (4.2) 

fc>0 

is bounded. In components, we calculate using the expression of ad^ that 

{Z l)u = V U ( f x k ~ X x i : (4-3) 
exp(z(A fc - A^)) - 1 
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Note that for any bounded operator A and -B, we always have 
||ad A (S)|| a < 2C a \\A\\jB\\ a 
where C a is given by Lemma 4.1 We define now the following numbers: 

Co = ft and Ci = 2C a \\Zt\\ , for £ > 1. 
Using (3.2) and Lemma 4.3, we easily see that we have the estimates 

v/>i, 2^-(/+i)Ofi<ii^ii a E E ^■•■C4- 

a fc>o ' 4+-+4=^ 

Now for any p such that tt < p < 2tt, there exist a constant M such that for all 
k,\B k \< k\Mp~ k . Hence we can write 

V*>1, ^-(£+l)0 +1 <Af||F|| a X; E C4---C4- 

a fc>0 £i+---+£ k =£ 

Let be the formal series £(£) = Y^£>o^ Ci- Multiplying the previous equation 
by t and summing over £ > 0, we find 

^r?(t)<M\\V\\ Y P ~ k ((t) k = M\\V\\ - — Ij-^. 
2C a s w " "« 1-C(t)/P 

Let r/(i) be the solution of the differential equation: 

r,'(t) = 2MC a \\V\\ a l _^ w „(0)=*. 
Taking p = 37r/2, we easily see that for t < 32 mc \\v\\ ' ^ e som tion can be 

Oi 

written 



and defines an analytic function of t. Expanding rj(t) = X^x)*^' we see that 
the coefficients satisfy the relations rjo = vr and 

V£>1, ^r(^ + l)%+i= E ^' E 

a fc>o fi+-+4=^ 

with p=^-. By induction, this shows that ^ < 77^. Moreover, for all z S C with 
|z| < 32 mc \\v\\ ' we nave as the coefficients are positive, 



IC(*)I 



E^ 



£=0 



<E^ = c(W)<i7(W)<^- 
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Using Cauchy estimates, we see that 
V€>1, \\Ze\\ 



1 ^ 1 C W (0) 3tt /32MC a ||y|| ^ 
Q = — — < 



2C a ~ 2C a l\ ~ 4C, 
The theorem is now proved by setting 

V(h) = Zi, and = //~ 2 ^ 



7T 



^>2 



32MC Q 



The estimate 



which defines a convergent power series for \h\ < ho 

(2.1) on V(h) is then an easy consequence of (4.2). The estimate (2.1) on W(h) 
is easily proved. 



5 Modified energy 

We give now the proof of Proposition 2.3. 
For all x G K, we have 



r y 2 

arctan(x) — x = — ^dy. 

Jo i + ir 



For A; G Z d , this yields 



- arctan (- L - L ) - fc 2 = / - n dy. 

h K 2 > 1 1 hj l + y 2 y 



Let 7 G [0, 2] , it is clear that for all y G 

,2 



l + y 2 



< u 7 . 



Hence we have for all k G Z Q 



2 + ^1*1% IH2 

— arctan ( — - — J — 



This shows that for all v, 



2 r^|fc| 2 /2 
< -r / y 7 dy < C7i 7 |Af 7+2 . 
h Jo 



(f | — — arctan {—^-)\v) — \v\ — A\v) 



< Ch?\ 



I #1+7 ■ 



Now we have 



(v | V(h) | v) - (v | V\v) = ^2^-{v I »*a4 o(fc) (V) 



fe>i 



(5.1) 
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Recall that Zo(h) = — 2arctan (^-) is a positive operator. The operator Zo(Zi) 1 / 2 
is hence well defined, and for an operator W we have in components 

(Z„(/0 1/2 WOm = (2arctan )) V Vfc<. 
Hence we have for all a > 1, 

||Z (/0 1/2 Wi Q < v^H^IL and ||WZ W 1/2 || Q < \^||W|| Q • 

Now using Lemma 4.2 and the fact that Zq(K) is symmetric, we have for all v 
and all operator W 

\(v\ a d Zo{h) (W)\v)\ < (llZoih^W^+llWZoih^JlZoih) 1 ^ \\v\\ 

< 2>/7r||W|| a \\Z {h) 1/2 v\\ \\v\\ . 

Hence we have 



\(v | V(h) \v)-(v\V\v)\<2j2 ^7r fc_1/2 ||^|| Q ||Z (/i) 1/2 w|| ||u| 



k>l 



< C\\V\\ a ||Z (/i) 1/2 ?;|| 

Using (5.1) with 7 = 0, this shows that 

\(v\V{h)\v) - (v\V\v) \ < C\\V\\ a h\\u\\ H1 \\u\\ . 
Finally, we easily have using (2.1) that 

\(v\W{h)\v)\ < C\\V\\ a h\\u\\ 2 . 
Summing the previous inequalities with 7 = (5 in (5.1) we have that 

(u\S(h)\u) -(u\A + V\u) < Ch^\\uf Hl+0 + C\\V\\ a h\\u\\ Hl ||u|| 
and this yields the result. 

6 Bounds for the numerical solution 

We prove now Corollary 2.4. Note that Eqn. (2.4) is classic. 

Using the fact that V is symmetric, we have for all n, ||ti n || = ||it°|| where 
II • || denotes the L 2 norm. 
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Using Lemma 4.2, we can write for all v £ L 2 , 

{v\S(h)\v) = t{v\-2 arctan {^—) \v) + (v\V{h) + hW(h) \ v) 
whence using (2.1), Lemma 4.2 and the fact that Z is a positive operator 

IH S(h)\v)\ > ±(v\ - 2arctan (^>> - C\\V\\ a \\vf . 
Hence using (2.4) we have that for all n, 

\{u n \ -2 arctan {\-)\u n ) < (u n \S(h)\u n ) + C\\V\\ \\u n f 

lb _ 

< (^isc/i)!^ ) + CH^H^ ||^°|| 2 . 

Using (2.3) with f3 = 0, we find that there exists a constant such that for all n, 
Uu n \ -2 arctan (— )\u n ) < C \\u f Hl . (6.1) 

Now we have for all x > 

1 ,lv , 1 2x , . 
x > — ==>- arctan x > arctan ( — ) and x < — =^ arctan x > — . (6.2) 

2 y 2 J ~ 2 3 v ; 

Applying this inequality to (6.1) by considering the set of frequencies h\k\ 2 < 1 
and h\k\ 2 > 1 immediately yields the result. 



7 Higher order approximations 

In this section we further investigate the long time behaviour by numerical sim- 
ulations and consider higher-order numerical schemes. 

We perform the simulations with d = 1, u° = 2/(2 — cos(x)) and V(x) = 
cos(x) + sin(6x). In the next figures, we show the maximal size of the oscillations 
of the truncated H 1 norm 

20 1/2 

( ^ (1 + |Af)K| 2 ) (7.1) 

k=-20 

along the numerical solution u n from t = to t = 50, and for stepsize ranging 
from h = 0.01 to h = 0.1. 

As expected, we see that this quantity is uniformly bounded for the splitting 
scheme (1.5) (Figure 1). 
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Figure 1: Midpoint approximation of the exponential. 



As explained in Remark 2.6, our methods easily extends to the Strang splitting 
scheme (2.6). Considering the alternative Strang splitting 

R{-ihA/2) exp(-ihV)R(-ihA/2), 

the same argument does not apply straightforwardly. The obstruction occurs in 
Lemma 4.3 where R(—ihA) is replaced by R(—ihA/2) 2 in the definition of the 
operator Zq, transforming ir by 2ir in inequality (4.1). 

Nevertheless, as shown in Figure 2, the same uniform conservation phenomenon 
can be observed. This might be justified using the fact that the operator Z\ de- 
fined in (4.3) still makes sense in this situation. 

Next we consider schemes of the form 

s 

exp{ihV) Yl R{-ljhA) (7.2) 
i=i 

where jj & R, j = 1, . . . , s are coefficients satisfying 71 + . ■ ■ + 7 S = 1. Such an 
approximation will be a higher order approximation of the splitting scheme (1.3) 
for suitable jj satisfying given algebraic conditions (see for instance [8, Chap 
III]). Of course, all these schemes remain symplectic and preserve the L 2 norm. 




0.04 



c > < > 3 5 



Figure 2: Strang splitting R(-ihA/2) exp(-ihV)R(-ihA/2). 
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In Figures 3, 4 and 5, we consider successively classical symmetric composition 
methods of order 4, 6 and 8 (see [8, Chap V] and the references therein). The 
method of order 4 is the triple jump method for which s = 3, 

1 2 x /3 
71=73 = ^^, and l2 = -—^. (7.3) 

The methods of order 6 corresponds to the methods given by Yoshida (see [15] 
and [8, Section V.3.2]) and requires s = 7, while the method of order 8 is the 
methods given by Suzuki & Umeno, see [14], and requires s = 15. 

What we observe is that for the method of order 4, the situation is similar 
to the previous cases (regularity conservation), but for the methods of order 6 
and 8, resonances appear: for specific values of the stepsize, the regularity of the 
numerical solution deteriorates. 




0.06 



0.08 



Figure 3: Order 4 approximation of the exponential. 



JL, 



C > « > 
C > ? 5 

c > < > 

C > - 1 

0.3 



0.04 



Figure 4: Order 6 approximation of the exponential. 

Finally, we plot in Figure 6 the same simulation for the "exact" splitting 
scheme (1.5). In this last situation, it is known that the resonances appear for 
step sizes h such that h(k 2 — I 2 ) is close to a multiple of 2n for some k and I € Z 
(see [3]). 

The fact that the method of order 4 possesses a modified energy can easily 
seen: With the values of 71, 72 and 73 given in (7.3), we have 

R(-j 1 hA)R(--y 2 hA)R(--f 3 hA) = exp(iZ ) 
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c > < > ? s 



Figure 5: Order 8 approximation of the exponential. 
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Figure 6: Exact splitting. 

where 

Z = -4 arctan (^^gj) + 2 arctan (^^y) = -G(*A/2) (7.4) 
with 

G{x) = 4 arctan ( 2 _ 2l/3 ) - 2 arctan ( 2 _ 2l/3 )- 

It is easy to see that for all a; > G(x) is an increasing function such that 
G(x) G (0, 7r). Hence Lemma 4.3 remains valid for this Zq. Using the same 
techniques as before, and bounds like (6.2) still valid for the function G(x), we 
can show the existence of a modified energy for this method, explaining the 
absence of resonances. 

Note that in the same spirit, we could consider symmetric composition meth- 
ods based on the order two Strang splitting (2.6) to build higher order methods 
of the form 

s 

J exp(i-/ j hV/2)R(-i~f j hA) exp(ijjhV/2) 

3=1 

to approximate (1.1). A general strategy to show the existence of a modified 
energy for this method would be to search for an operator Z(t) such that for all 
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t > 0, 

s 

exp(iZ(t)) = Y[e*p(hjtV/2)R(-i~f j hA)exp(i~f j tV/2) 

with 

s 

Zq = — 2 arctan(/t7jA/2). 

In the case of the triple jump method, this operator can be written (7.4), and the 
same argument as above shows the existence of a modified energy for this method 
by using the same kind of techniques. We do not give the details here. The 
derivation of higher order methods possessing a modified energy is an interesting 
question that will be addressed in future studies. 
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